%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This adapted code comes from "Identification Properties for 
% Estimating the Impact of Regulation on Markups and Productivity".
% by Sampi, James; Jooste, Charl; Vostroknutova, Ekaterina, January 2021.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Table_C1_b_c (pooled and yearly effects)

% Note that the code `markups_sampi_estimation.m` is a single code that appends (in order) 
% all MATLAB codes used. This sequence uses three initial inputs: 
% `like.m` and `standardizeT_new.m` MATLAB functions from Sampi et al, and `analysis_data.dta` 
% (events data: new suppliers to multinationals). It produces all parts of 
% Supplementary_Table_C1_b_c (pooled and yearly effects).

% `markups_sampi_estimation_global.m`: implements Sampi's et al. (2021) original method.
% It produces one markup impact coefficient after event and its standard deviation. 
% It produces the output: `Supplementary_Table_C1_b_pooled.csv`. 

% `markups_sampi_estimation_global_never.m`: implements Sampi's et al. (2021) method, however,
% it excludes never-suppliers as controls and uses only eventually-treated. It produces the 
% output: `Supplementary_Table_C1_c_pooled.csv`. 

% `markups_sampi_estimation_rel_med_aleat.m`: implements Sampi's et al. (2021) method modified
% to separate the markup impact coefficient by event-relative time horizon (years before and 
% after event). It produces eight coefficients (one by relative horizon and excluding relative
% time $t=-1$ since it is used as the reference base) and their standard deviation. It produces
% the output: `Supplementary_Table_C1_b_yearly.csv`. 

% `markups_sampi_estimation_rel_med_never_aleat.m`: implements Sampi's et al. (2021) method 
% modified to separate the markup impact coefficient by event-relative time horizon. However, 
% it excludes never-suppliers as controls and uses only eventually-treated. It produces the 
% output: `Supplementary_Table_C1_c_yearly.csv`. 


%%
%markups_sampi_estimation_global (Table_C1_b_pooled)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear  
close all
clc 

rng (1)

%adding folders to path
cd ..\.. %only run in first try
addpath 'codes\all_codes\'

Kr = readmatrix('temp/matlab_input_capital_global.csv');    %capital
L = readmatrix('temp/matlab_input_trabaj_global.csv');      %workers
s = readmatrix('temp/matlab_input_labor_sh_2_global.csv'); %labor share
Yr = readmatrix('temp/matlab_input_ventas_global.csv');     %sales
Dummy = readmatrix('temp/matlab_input_D_global.csv');       %event dummy

Kr(Kr()==1.123456) = NaN;
L(L()==1.123456) = NaN;
s(s()==1.123456) = NaN;
Yr(Yr()==1.123456) = NaN;

tic
     
N = size(L,1); %number of firms
T = size(L,2); %time dimension

cobb = 1; % 1 for Cobb-Douglas production function, 0 otherwise - CES 

% Model estimation %

y1 = -log(s); %the negative of the log (labor share)
y2 = log(Yr); % the log of total sales 

x2 = log(Kr); % log of nominal stock of capital
x3 = log(L); % log of employees 
x4 = log(log(Kr)-log(L)); % proxy for elas as stated in Equation 31 in the paper

% Estimation using Q-MLE 
options = optimset('Display','off','TolX',1e-6);
options.LargeScale='off';
options.Algorithm = 'bfgs';
options.MaxFunEvals = Inf;

theta0    = 0.1*ones(13,1);   %Initial parameters  
if cobb==1
    theta0(10)=0; %if CB then output elasticity is constant
else
    theta0(10)=1; % if CES then output elasticity changes over time as per Equation 31
end

[x] = ...
fminunc(@(theta0) like(theta0,y1,y2,x2,x3,x4,Dummy),theta0,options);

%main parameters of interest
coefMarkup = x(11);

%bootstrap
bootfun = @(y1,y2,x2,x3,x4,Dummy) fminunc(@(theta0) like(theta0,y1,y2,x2,x3,x4,Dummy),theta0,options);  
se = std(bootstrp(500,bootfun,y1,y2,x2,x3,x4,Dummy));

%saving coef vector
final_vec = [coefMarkup se(11)];
Folder = 'codes\results\3-appendix_c\';
File = fullfile(Folder, sprintf('Supplementary_Table_C1_b_pooled.csv'));
writematrix(final_vec, File);

toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%
%markups_sampi_estimation_global_never (Table_C1_c_pooled)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear  
close all
clc 

rng (1)

%adding folders to path
cd ..\.. %only run in first try
addpath 'codes\all_codes\'

Kr = readmatrix('temp/matlab_input_capital_global_never.csv');    %capital
L = readmatrix('temp/matlab_input_trabaj_global_never.csv');      %workers
s = readmatrix('temp/matlab_input_labor_sh_2_global_never.csv'); %labor share
Yr = readmatrix('temp/matlab_input_ventas_global_never.csv');     %sales
Dummy = readmatrix('temp/matlab_input_D_global_never.csv');       %event dummy

Kr(Kr()==1.123456) = NaN;
L(L()==1.123456) = NaN;
s(s()==1.123456) = NaN;
Yr(Yr()==1.123456) = NaN; 
Dummy(Dummy()==1.123456) = NaN; 

tic
     
N = size(L,1); %number of firms
T = size(L,2); %time dimension

cobb = 1; % 1 for Cobb-Douglas production function, 0 otherwise - CES 

% Model estimation %

y1 = -log(s); %the negative of the log (labor share)
y2 = log(Yr); % the log of total sales 

x2 = log(Kr); % log of nominal stock of capital
x3 = log(L); % log of employees 
x4 = log(log(Kr)-log(L)); % proxy for elas as stated in Equation 31 in the paper 

% Estimation using Q-MLE 
options = optimset('Display','off','TolX',1e-6);
options.LargeScale='off';
options.Algorithm = 'bfgs';
options.MaxFunEvals = Inf;

theta0    = 0.1*ones(13,1);   %Initial parameters  
if cobb==1
    theta0(10)=0; %if CB then output elasticity is constant
else
    theta0(10)=1; % if CES then output elasticity changes over time as per Equation 31
end

[x] = ...
fminunc(@(theta0) like(theta0,y1,y2,x2,x3,x4,Dummy),theta0,options);

%main parameters of interest
coefMarkup = x(11);

%bootstrap
bootfun = @(y1,y2,x2,x3,x4,Dummy) fminunc(@(theta0) like(theta0,y1,y2,x2,x3,x4,Dummy),theta0,options);  
se = std(bootstrp(500,bootfun,y1,y2,x2,x3,x4,Dummy));

%saving coef vector
final_vec = [coefMarkup se(11)];
Folder = 'codes\results\3-appendix_c\';
File = fullfile(Folder, sprintf('Supplementary_Table_C1_c_pooled.csv'));
writematrix(final_vec, File);

toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%
%markups_sampi_estimation_rel_med_aleat (Table_C1_b_yearly)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear  
close all
clc 

rng (1)

%adding folders to path
cd ..\.. %only run in first try
addpath 'codes\all_codes\'

Kr_base = readmatrix('temp/matlab_input_capital_rel_med_aleat.csv');    %capital
L_base = readmatrix('temp/matlab_input_trabaj_rel_med_aleat.csv');      %workers
s_base = readmatrix('temp/matlab_input_labor_sh_2_rel_med_aleat.csv'); %labor share
Yr_base = readmatrix('temp/matlab_input_ventas_rel_med_aleat.csv');     %sales
Dummy_base = readmatrix('temp/matlab_input_D_rel_med_aleat.csv');       %event dummy

Kr_base(Kr_base()==1.123456) = NaN;
L_base(L_base()==1.123456) = NaN;
s_base(s_base()==1.123456) = NaN;
Yr_base(Yr_base()==1.123456) = NaN; 
Dummy_base(Dummy_base()==1.123456) = NaN; 

final_vec = [];
N = size(L_base,1); %number of firms

tic

for t = [1, 2, 3, 5, 6, 7, 8, 9]   
           
if t<=3       
  Kr = Kr_base(:, (t):(4));
  L = L_base(:, (t):(4));
  s = s_base(:, (t):(4));
  Yr = Yr_base(:, (t):(4));
  Dummy = Dummy_base(:, [t:3 5]);
else
  Kr = Kr_base(:, (4):(t));
  L = L_base(:, (4):(t));
  s = s_base(:, (4):(t));
  Yr = Yr_base(:, (4):(t));
  Dummy = Dummy_base(:, (4):(t));
end


T = size(L,2); %time dimension

cobb = 1; % 1 for Cobb-Douglas production function, 0 otherwise - CES 

% Model estimation %

y1 = -log(s); %the negative of the log (labor share)
y2 = log(Yr); % the log of total sales 

x2 = log(Kr); % log of nominal stock of capital
x3 = log(L); % log of employees 
x4 = log(log(Kr)-log(L)); % proxy for elas as stated in Equation 31 in the paper 

% Estimation using Q-MLE 
options = optimset('Display','off','TolX',1e-6);
options.LargeScale='off';
options.Algorithm = 'bfgs';
options.MaxFunEvals = Inf;

theta0    = 0.1*ones(13,1);   %Initial parameters  
if cobb==1
    theta0(10)=0; %if CB then output elasticity is constant
else
    theta0(10)=1; % if CES then output elasticity changes over time as per Equation 31
end

[x] = ...
fminunc(@(theta0) like(theta0,y1,y2,x2,x3,x4,Dummy),theta0,options);

%main parameters of interest
coefMarkup = x(11);

%bootstrap
bootfun = @(y1,y2,x2,x3,x4,Dummy) fminunc(@(theta0) like(theta0,y1,y2,x2,x3,x4,Dummy),theta0,options);  
se = std(bootstrp(500,bootfun,y1,y2,x2,x3,x4,Dummy));

%saving coef vector
vec1 = [ (t-5) coefMarkup se(11)];
final_vec = [final_vec; vec1];

end
 
Folder = 'codes\results\3-appendix_c\';
File = fullfile(Folder, sprintf('Supplementary_Table_C1_b_yearly.csv'));
writematrix(final_vec, File);
 
toc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%
%markups_sampi_estimation_rel_med_never_aleat (Table_C1_c_yearly)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear  
close all
clc 

rng (1)

%adding folders to path
cd ..\.. %only run in first try
addpath 'codes\all_codes\'

Kr_base = readmatrix('temp/matlab_input_capital_rel_med_never_aleat.csv');    %capital
L_base = readmatrix('temp/matlab_input_trabaj_rel_med_never_aleat.csv');      %workers
s_base = readmatrix('temp/matlab_input_labor_sh_2_rel_med_never_aleat.csv'); %labor share
Yr_base = readmatrix('temp/matlab_input_ventas_rel_med_never_aleat.csv');     %sales
Dummy_base = readmatrix('temp/matlab_input_D_rel_med_never_aleat.csv');       %event dummy

Kr_base(Kr_base()==1.123456) = NaN;
L_base(L_base()==1.123456) = NaN;
s_base(s_base()==1.123456) = NaN;
Yr_base(Yr_base()==1.123456) = NaN; 
Dummy_base(Dummy_base()==1.123456) = NaN; 

final_vec = [];
N = size(L_base,1); %number of firms

tic

for t = [1, 2, 3, 5, 6, 7, 8, 9]   
           
if t<=3       
  Kr = Kr_base(:, (t):(4));
  L = L_base(:, (t):(4));
  s = s_base(:, (t):(4));
  Yr = Yr_base(:, (t):(4));
  Dummy = Dummy_base(:, [t:3 5]);
else
  Kr = Kr_base(:, (4):(t));
  L = L_base(:, (4):(t));
  s = s_base(:, (4):(t));
  Yr = Yr_base(:, (4):(t));
  Dummy = Dummy_base(:, (4):(t));
end


T = size(L,2); %time dimension

cobb = 1; % 1 for Cobb-Douglas production function, 0 otherwise - CES 

% Model estimation %

y1 = -log(s); %the negative of the log (labor share)
y2 = log(Yr); % the log of total sales 

x2 = log(Kr); % log of nominal stock of capital
x3 = log(L); % log of employees 
x4 = log(log(Kr)-log(L)); % proxy for elas as stated in Equation 31 in the paper 

% Estimation using Q-MLE 
options = optimset('Display','off','TolX',1e-6);
options.LargeScale='off';
options.Algorithm = 'bfgs';
options.MaxFunEvals = Inf;

theta0    = 0.1*ones(13,1);   %Initial parameters  
if cobb==1
    theta0(10)=0; %if CB then output elasticity is constant
else
    theta0(10)=1; % if CES then output elasticity changes over time as per Equation 31
end

[x] = ...
fminunc(@(theta0) like(theta0,y1,y2,x2,x3,x4,Dummy),theta0,options);

%main parameters of interest
coefMarkup = x(11);

%bootstrap
bootfun = @(y1,y2,x2,x3,x4,Dummy) fminunc(@(theta0) like(theta0,y1,y2,x2,x3,x4,Dummy),theta0,options);  
se = std(bootstrp(500,bootfun,y1,y2,x2,x3,x4,Dummy));

%saving coef vector
vec1 = [ (t-5) coefMarkup se(11)];
final_vec = [final_vec; vec1];

end
 
Folder = 'codes\results\3-appendix_c\';
File = fullfile(Folder, sprintf('Supplementary_Table_C1_c_yearly.csv'));
writematrix(final_vec, File);
 
toc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%
% Matlab functions from Sampi et al (2021)

function y = standardizeT_new(x)

% Function to make your data have mean 0 and variance 1.
% Data in x are Txp, i.e. T time series observations times p variables
%y = (x - repmat(mean(x,1),size(x,1),1)); %cross-seciton mean
%y = (x - repmat(mean(x,2),1,size(x,2))); %time dimension mean
y = mean(mean(x, 'omitnan'), 'omitnan') + (x - repmat(mean(x,1, 'omitnan'),size(x,1),1)-repmat(mean(x,2, 'omitnan'),1,size(x,2)));
end

function LL=like(x,y1,y2,x2,x3,x4,D)    

    y1 = standardizeT_new(y1);
    y2 = standardizeT_new(y2);
    x2 = standardizeT_new(x2);
    x3 = standardizeT_new(x3);
    x4 = standardizeT_new(x4);

    D = standardizeT_new(D);

    [T, N]= size(y1');
    
    Q1 = exp(log(x(1)));
    FF = inv(Q1);
    
    K = x2; 
    Y = y2; 
    
    for t = 2:T        
        nonpara(:,t) =  x(3)*K(:,t) + x(4)*x3(:,t) + ...
            + x(5)*K(:,t).*x3(:,t) + ...
            + 0.5*x(6)*K(:,t).^2 + 0.5*x(7)*x3(:,t).^2;       
        
        e1(:,t) = y1(:,t)  - x(10)*x4(:,t) - x(11)*D(:,t); 
        e2(:,t) = Y(:,t) - nonpara(:,t) - x(12)*(Y(:,t-1) - nonpara(:,t-1)) - x(13)*D(:,t);
        
        e(:,t) = [e1(:,t); e2(:,t)];
        e(isnan(e)) = 0;
    %%%%%%%%%%%%%%%%%%%%% computed the joint probability of the log-likelihood function
        ly(t) = sum(e(:,t)'*FF*e(:,t),2, 'omitnan');
        logl(t) = -0.5*ly(t) -(2/2)*N*log(Q1);       
    end        
    
LL=-sum(logl(2:end), 'omitnan');
end
